function savings_EWM_onedim = generate_results_fixedshare_onedimension(bsimax,sample,share,fix_range,wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase)
% This function calculates cost savings or energy reduction from EWM
% onedimensional rules under a fixed budget using grid search. 

% Inputs: 
% (1) bsimax: number of bootstrap runs
% (2) sample: cross-sectional dataset with electricity consumption, 
% covariates, and propensity scores
% (3) share: share of households allowed to treat under the fixed budget
% (4) fix_range: percentage point variation from the share that is allowed
% (5) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0)
% (6) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (7) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (8) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (9) SMC: =1 use social marginal cost; =0 use retail electricity price
% (10) savedir: directory to save the savings table 
% (11) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) savings_EWM_onedim: a table summarizing the savings point estimate,
% share of households to treat

%% Gather inputs

[X,Y,D,n,ps]=generate_inputs_quadrant_rule(sample,wave,tcost,covariate,SMC);

%% Estimate g using IPW
g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

%% Level search
% sort rows of X -> [Xr, gr]
[Xr, Ind] = sortrows(X);
gr = g(Ind);

% descretize baseline usage
dx1=10;
[N,boo1,xx1]=histcounts(Xr(:,2),...
    'BinWidth',[dx1],'Normalization','count'); %

nu = length(boo1);

Xu=nan(nu,1); gu=nan(nu,1); nw = nan(nu,1);
k=1;
for i = 1:length(boo1)
    bin_boo = logical((xx1==i));
    % identify individuals in each level of baseline consumption
    Xu(k,:)=boo1(i);
    gu(k,1) = sum(gr(bin_boo)); % net savings for each unique covariate combo
    nw(k,1) = sum(bin_boo); % count number of individuals in the pair of covariate combo
    k=k+1;
end

x1 = Xu(:,1);

% Create grids of midpoints between distinct values of x1/x2
u_x1 = boo1';

eps_x1 = min(u_x1(2:end) - u_x1(1:end-1))/2;
grid_x1 = [u_x1(1)-eps_x1; (u_x1(2:end) + u_x1(1:end-1))/2; u_x1(end)+eps_x1];
k1 = size(grid_x1,1);

%%% hardcoded for 1D
k2=1; Ndim=1;

% Grid search
tic1= tic;
c_share_candidate=cell(k1,k2);
myfun_okayShare = @(x) x-share>-fix_range; % define the fix share constraint
fractn=nan(k1,1);
m_W = inf(k1,k2,Ndim*2);
minW = min(m_W);

for i1=1:k1
    boo_i1=(x1 - grid_x1(i1));
    temp_ratio = zeros(1,Ndim*2);
    for sign1 = -1:2:1 %for i = [-1,1]
        select = ((sign1.*boo_i1)>=0);
        treat_ratio_candidate = sum(nw.*select)/n;
        if sign1==-1
            idx_sign = 1; % get the index to save the signs 
        else
            idx_sign = 2;
        end
        temp_ratio(idx_sign) = treat_ratio_candidate;
        %
        if myfun_okayShare(treat_ratio_candidate) % only calculate W for the rules that meet the fix share constraint
            fraction = min(1,share./treat_ratio_candidate);
            W_rationed = sum(gu.*select).*fraction;
            m_W(i1,1,idx_sign) = W_rationed;
            fractn (i1,1) = fraction;
        else
            W_rationed = inf; % otherwise W remains to be inf
        end
        
        if (W_rationed < minW)
            min_i1 = i1;
            min_sign1 = sign1;
            minW = W_rationed;
        end
    end
    c_share_candidate{i1,1} = temp_ratio;
end

candidate_grid=cellfun(myfun_okayShare,c_share_candidate,'UniformOutput',0);
N_Rule_okayShare = sum([candidate_grid{:}],'all'); % number of rules that meet the fix share constraint

fprintf('Serach takes %.2f sec\n',toc(tic1))

%% Save coefficients and other variables for graphs

in_Ghat = (min_sign1*(x1 - grid_x1(min_i1))>=0);
select_fraction=fractn(min_i1,1);
percent=sum(nw.*in_Ghat)/n*select_fraction;
savings=sum(gu.*in_Ghat)/n*select_fraction;

switch winter_prevuse 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 % cost_savings
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end
s_range = sprintf('fix%0.2f',fix_range*100);
filename_coefs=sprintf('coef_fixedshare_onedim_%s_%s%s_wave%1.0f_%s.mat',covariate,s_cost,s_winter,wave,s_range);
filename = sprintf('%s_%s',fbase,filename_coefs);
fpathf = fullfile(savedir,filename);
fprintf('Saving "%s" to "%s" ... ',filename,savedir); % must use newline symbol "\n" at the end
save(fpathf,'min_i1','min_sign1',...
    'in_Ghat','Xu','nw','gu','n','covariate','minW',...
    'percent','savings','c_share_candidate','candidate_grid','N_Rule_okayShare','fix_range','share','fractn','select_fraction');

%% Output for table
savings_EWM_onedim=nan(1,4);
savings_EWM_onedim=array2table(savings_EWM_onedim,'VariableNames',{'rules','Covariates','Share\ treated','Net\ cost\ changes'});
format shortg
savings_EWM_onedim.(1)=  {'quadrant'};
savings_EWM_onedim.(2) = {covariate};
savings_EWM_onedim.(3) = round(percent*100,2);
savings_EWM_onedim.(4) = round(savings,3);
end